home *** CD-ROM | disk | FTP | other *** search
/ OpenGL Superbible (2nd Edition) / OpenGL SuperBible e2.iso / tools / GLUT-3.7 / LIB / GLE / VVECTOR.H < prev   
Encoding:
C/C++ Source or Header  |  1998-08-12  |  34.2 KB  |  1,340 lines
  1.  
  2. /*
  3.  * vvector.h
  4.  *
  5.  * FUNCTION:
  6.  * This file contains a number of utilities useful for handling
  7.  * 3D vectors
  8.  * 
  9.  * HISTORY:
  10.  * Written by Linas Vepstas, August 1991
  11.  * Added 2D code, March 1993
  12.  * Added Outer products, C++ proofed, Linas Vepstas October 1993
  13.  */
  14.  
  15. #ifndef __GUTIL_VECTOR_H__
  16. #define __GUTIL_VECTOR_H__
  17.  
  18. #if defined(__cplusplus) || defined(c_plusplus)
  19. extern "C" {
  20. #endif
  21.  
  22.  
  23. #include <math.h>
  24. #include "port.h"
  25.  
  26. /* ========================================================== */
  27. /* Zero out a 2D vector */
  28.  
  29. #define VEC_ZERO_2(a)                \
  30. {                        \
  31.    (a)[0] = (a)[1] = 0.0;            \
  32. }
  33.  
  34. /* ========================================================== */
  35. /* Zero out a 3D vector */
  36.  
  37. #define VEC_ZERO(a)                \
  38. {                        \
  39.    (a)[0] = (a)[1] = (a)[2] = 0.0;        \
  40. }
  41.  
  42. /* ========================================================== */
  43. /* Zero out a 4D vector */
  44.  
  45. #define VEC_ZERO_4(a)                \
  46. {                        \
  47.    (a)[0] = (a)[1] = (a)[2] = (a)[3] = 0.0;    \
  48. }
  49.  
  50. /* ========================================================== */
  51. /* Vector copy */
  52.  
  53. #define VEC_COPY_2(b,a)                \
  54. {                        \
  55.    (b)[0] = (a)[0];                \
  56.    (b)[1] = (a)[1];                \
  57. }
  58.  
  59. /* ========================================================== */
  60. /* Copy 3D vector */
  61.  
  62. #define VEC_COPY(b,a)                \
  63. {                        \
  64.    (b)[0] = (a)[0];                \
  65.    (b)[1] = (a)[1];                \
  66.    (b)[2] = (a)[2];                \
  67. }
  68.  
  69. /* ========================================================== */
  70. /* Copy 4D vector */
  71.  
  72. #define VEC_COPY_4(b,a)                \
  73. {                        \
  74.    (b)[0] = (a)[0];                \
  75.    (b)[1] = (a)[1];                \
  76.    (b)[2] = (a)[2];                \
  77.    (b)[3] = (a)[3];                \
  78. }
  79.  
  80. /* ========================================================== */
  81. /* Vector difference */
  82.  
  83. #define VEC_DIFF_2(v21,v2,v1)            \
  84. {                        \
  85.    (v21)[0] = (v2)[0] - (v1)[0];        \
  86.    (v21)[1] = (v2)[1] - (v1)[1];        \
  87. }
  88.  
  89. /* ========================================================== */
  90. /* Vector difference */
  91.  
  92. #define VEC_DIFF(v21,v2,v1)            \
  93. {                        \
  94.    (v21)[0] = (v2)[0] - (v1)[0];        \
  95.    (v21)[1] = (v2)[1] - (v1)[1];        \
  96.    (v21)[2] = (v2)[2] - (v1)[2];        \
  97. }
  98.  
  99. /* ========================================================== */
  100. /* Vector difference */
  101.  
  102. #define VEC_DIFF_4(v21,v2,v1)            \
  103. {                        \
  104.    (v21)[0] = (v2)[0] - (v1)[0];        \
  105.    (v21)[1] = (v2)[1] - (v1)[1];        \
  106.    (v21)[2] = (v2)[2] - (v1)[2];        \
  107.    (v21)[3] = (v2)[3] - (v1)[3];        \
  108. }
  109.  
  110. /* ========================================================== */
  111. /* Vector sum */
  112.  
  113. #define VEC_SUM_2(v21,v2,v1)            \
  114. {                        \
  115.    (v21)[0] = (v2)[0] + (v1)[0];        \
  116.    (v21)[1] = (v2)[1] + (v1)[1];        \
  117. }
  118.  
  119. /* ========================================================== */
  120. /* Vector sum */
  121.  
  122. #define VEC_SUM(v21,v2,v1)            \
  123. {                        \
  124.    (v21)[0] = (v2)[0] + (v1)[0];        \
  125.    (v21)[1] = (v2)[1] + (v1)[1];        \
  126.    (v21)[2] = (v2)[2] + (v1)[2];        \
  127. }
  128.  
  129. /* ========================================================== */
  130. /* Vector sum */
  131.  
  132. #define VEC_SUM_4(v21,v2,v1)            \
  133. {                        \
  134.    (v21)[0] = (v2)[0] + (v1)[0];        \
  135.    (v21)[1] = (v2)[1] + (v1)[1];        \
  136.    (v21)[2] = (v2)[2] + (v1)[2];        \
  137.    (v21)[3] = (v2)[3] + (v1)[3];        \
  138. }
  139.  
  140. /* ========================================================== */
  141. /* scalar times vector */
  142.  
  143. #define VEC_SCALE_2(c,a,b)            \
  144. {                        \
  145.    (c)[0] = (a)*(b)[0];                \
  146.    (c)[1] = (a)*(b)[1];                \
  147. }
  148.  
  149. /* ========================================================== */
  150. /* scalar times vector */
  151.  
  152. #define VEC_SCALE(c,a,b)            \
  153. {                        \
  154.    (c)[0] = (a)*(b)[0];                \
  155.    (c)[1] = (a)*(b)[1];                \
  156.    (c)[2] = (a)*(b)[2];                \
  157. }
  158.  
  159. /* ========================================================== */
  160. /* scalar times vector */
  161.  
  162. #define VEC_SCALE_4(c,a,b)            \
  163. {                        \
  164.    (c)[0] = (a)*(b)[0];                \
  165.    (c)[1] = (a)*(b)[1];                \
  166.    (c)[2] = (a)*(b)[2];                \
  167.    (c)[3] = (a)*(b)[3];                \
  168. }
  169.  
  170. /* ========================================================== */
  171. /* accumulate scaled vector */
  172.  
  173. #define VEC_ACCUM_2(c,a,b)            \
  174. {                        \
  175.    (c)[0] += (a)*(b)[0];            \
  176.    (c)[1] += (a)*(b)[1];            \
  177. }
  178.  
  179. /* ========================================================== */
  180. /* accumulate scaled vector */
  181.  
  182. #define VEC_ACCUM(c,a,b)            \
  183. {                        \
  184.    (c)[0] += (a)*(b)[0];            \
  185.    (c)[1] += (a)*(b)[1];            \
  186.    (c)[2] += (a)*(b)[2];            \
  187. }
  188.  
  189. /* ========================================================== */
  190. /* accumulate scaled vector */
  191.  
  192. #define VEC_ACCUM_4(c,a,b)            \
  193. {                        \
  194.    (c)[0] += (a)*(b)[0];            \
  195.    (c)[1] += (a)*(b)[1];            \
  196.    (c)[2] += (a)*(b)[2];            \
  197.    (c)[3] += (a)*(b)[3];            \
  198. }
  199.  
  200. /* ========================================================== */
  201. /* Vector dot product */
  202.  
  203. #define VEC_DOT_PRODUCT_2(c,a,b)            \
  204. {                            \
  205.    c = (a)[0]*(b)[0] + (a)[1]*(b)[1];            \
  206. }
  207.  
  208. /* ========================================================== */
  209. /* Vector dot product */
  210.  
  211. #define VEC_DOT_PRODUCT(c,a,b)                \
  212. {                            \
  213.    c = (a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2]*(b)[2];    \
  214. }
  215.  
  216. /* ========================================================== */
  217. /* Vector dot product */
  218.  
  219. #define VEC_DOT_PRODUCT_4(c,a,b)            \
  220. {                            \
  221.    c = (a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2]*(b)[2] + (a)[3]*(b)[3] ;    \
  222. }
  223.  
  224. /* ========================================================== */
  225. /* vector impact parameter (squared) */
  226.  
  227. #define VEC_IMPACT_SQ(bsq,direction,position)        \
  228. {                            \
  229.    gleDouble len, llel;                    \
  230.    VEC_DOT_PRODUCT (len, position, position);        \
  231.    VEC_DOT_PRODUCT (llel, direction, position);        \
  232.    bsq = len - llel*llel;                \
  233. }
  234.  
  235. /* ========================================================== */
  236. /* vector impact parameter */
  237.  
  238. #define VEC_IMPACT(bsq,direction,position)        \
  239. {                            \
  240.    VEC_IMPACT_SQ(bsq,direction,position);        \
  241.    bsq = sqrt (bsq);                    \
  242. }
  243.  
  244. /* ========================================================== */
  245. /* Vector length */
  246.  
  247. #define VEC_LENGTH_2(len,a)            \
  248. {                        \
  249.    len = a[0]*a[0] + a[1]*a[1];            \
  250.    len = sqrt (len);                \
  251. }
  252.  
  253. /* ========================================================== */
  254. /* Vector length */
  255.  
  256. #define VEC_LENGTH(len,a)            \
  257. {                        \
  258.    len = (a)[0]*(a)[0] + (a)[1]*(a)[1];        \
  259.    len += (a)[2]*(a)[2];            \
  260.    len = sqrt (len);                \
  261. }
  262.  
  263. /* ========================================================== */
  264. /* Vector length */
  265.  
  266. #define VEC_LENGTH_4(len,a)            \
  267. {                        \
  268.    len = (a)[0]*(a)[0] + (a)[1]*(a)[1];        \
  269.    len += (a)[2]*(a)[2];            \
  270.    len += (a)[3] * (a)[3];            \
  271.    len = sqrt (len);                \
  272. }
  273.  
  274. /* ========================================================== */
  275. /* distance between two points */
  276.  
  277. #define VEC_DISTANCE(len,va,vb)            \
  278. {                        \
  279.     gleDouble tmp[4];                \
  280.     VEC_DIFF (tmp, vb, va);            \
  281.     VEC_LENGTH (len, tmp);            \
  282. }
  283.  
  284. /* ========================================================== */
  285. /* Vector length */
  286.  
  287. #define VEC_CONJUGATE_LENGTH(len,a)        \
  288. {                        \
  289.    len = 1.0 - a[0]*a[0] - a[1]*a[1] - a[2]*a[2];\
  290.    len = sqrt (len);                \
  291. }
  292.  
  293. /* ========================================================== */
  294. /* Vector length */
  295.  
  296. #define VEC_NORMALIZE(a)            \
  297. {                        \
  298.    double len;                    \
  299.    VEC_LENGTH (len,a);                \
  300.    if (len != 0.0) {                \
  301.       len = 1.0 / len;                \
  302.       a[0] *= len;                \
  303.       a[1] *= len;                \
  304.       a[2] *= len;                \
  305.    }                        \
  306. }
  307.  
  308. /* ========================================================== */
  309. /* Vector length */
  310.  
  311. #define VEC_RENORMALIZE(a,newlen)        \
  312. {                        \
  313.    double len;                    \
  314.    VEC_LENGTH (len,a);                \
  315.    if (len != 0.0) {                \
  316.       len = newlen / len;                \
  317.       a[0] *= len;                \
  318.       a[1] *= len;                \
  319.       a[2] *= len;                \
  320.    }                        \
  321. }
  322.  
  323. /* ========================================================== */
  324. /* 3D Vector cross product yeilding vector */
  325.  
  326. #define VEC_CROSS_PRODUCT(c,a,b)        \
  327. {                        \
  328.    c[0] = (a)[1] * (b)[2] - (a)[2] * (b)[1];    \
  329.    c[1] = (a)[2] * (b)[0] - (a)[0] * (b)[2];    \
  330.    c[2] = (a)[0] * (b)[1] - (a)[1] * (b)[0];    \
  331. }
  332.  
  333. /* ========================================================== */
  334. /* Vector perp -- assumes that n is of unit length 
  335.  * accepts vector v, subtracts out any component parallel to n */
  336.  
  337. #define VEC_PERP(vp,v,n)            \
  338. {                        \
  339.    double dot;                    \
  340.                         \
  341.    VEC_DOT_PRODUCT (dot, v, n);            \
  342.    vp[0] = (v)[0] - (dot) * (n)[0];        \
  343.    vp[1] = (v)[1] - (dot) * (n)[1];        \
  344.    vp[2] = (v)[2] - (dot) * (n)[2];        \
  345. }
  346.  
  347. /* ========================================================== */
  348. /* Vector parallel -- assumes that n is of unit length 
  349.  * accepts vector v, subtracts out any component perpendicular to n */
  350.  
  351. #define VEC_PARALLEL(vp,v,n)            \
  352. {                        \
  353.    double dot;                    \
  354.                         \
  355.    VEC_DOT_PRODUCT (dot, v, n);            \
  356.    vp[0] = (dot) * (n)[0];            \
  357.    vp[1] = (dot) * (n)[1];            \
  358.    vp[2] = (dot) * (n)[2];            \
  359. }
  360.  
  361. /* ========================================================== */
  362. /* Vector reflection -- assumes n is of unit length */
  363. /* Takes vector v, reflects it against reflector n, and returns vr */
  364.  
  365. #define VEC_REFLECT(vr,v,n)            \
  366. {                        \
  367.    double dot;                    \
  368.                         \
  369.    VEC_DOT_PRODUCT (dot, v, n);            \
  370.    vr[0] = (v)[0] - 2.0 * (dot) * (n)[0];    \
  371.    vr[1] = (v)[1] - 2.0 * (dot) * (n)[1];    \
  372.    vr[2] = (v)[2] - 2.0 * (dot) * (n)[2];    \
  373. }
  374.  
  375. /* ========================================================== */
  376. /* Vector blending */
  377. /* Takes two vectors a, b, blends them together */ 
  378.  
  379. #define VEC_BLEND(vr,sa,a,sb,b)            \
  380. {                        \
  381.                         \
  382.    vr[0] = (sa) * (a)[0] + (sb) * (b)[0];    \
  383.    vr[1] = (sa) * (a)[1] + (sb) * (b)[1];    \
  384.    vr[2] = (sa) * (a)[2] + (sb) * (b)[2];    \
  385. }
  386.  
  387. /* ========================================================== */
  388. /* Vector print */
  389.  
  390. #define VEC_PRINT_2(a)                    \
  391. {                            \
  392.    double len;                        \
  393.    VEC_LENGTH_2 (len, a);                    \
  394.    printf (" a is %f %f length of a is %f \n", a[0], a[1], len); \
  395. }
  396.  
  397. /* ========================================================== */
  398. /* Vector print */
  399.  
  400. #define VEC_PRINT(a)                    \
  401. {                            \
  402.    double len;                        \
  403.    VEC_LENGTH (len, (a));                \
  404.    printf (" a is %f %f %f length of a is %f \n", (a)[0], (a)[1], (a)[2], len); \
  405. }
  406.  
  407. /* ========================================================== */
  408. /* Vector print */
  409.  
  410. #define VEC_PRINT_4(a)                    \
  411. {                            \
  412.    double len;                        \
  413.    VEC_LENGTH_4 (len, (a));                \
  414.    printf (" a is %f %f %f %f length of a is %f \n",    \
  415.        (a)[0], (a)[1], (a)[2], (a)[3], len);        \
  416. }
  417.  
  418. /* ========================================================== */
  419. /* print matrix */
  420.  
  421. #define MAT_PRINT_4X4(mmm) {                \
  422.    int i,j;                        \
  423.    printf ("matrix mmm is \n");                \
  424.    if (mmm == NULL) {                    \
  425.       printf (" Null \n");                \
  426.    } else {                        \
  427.       for (i=0; i<4; i++) {                \
  428.          for (j=0; j<4; j++) {                \
  429.             printf ("%f ", mmm[i][j]);            \
  430.          }                        \
  431.          printf (" \n");                \
  432.       }                            \
  433.    }                            \
  434. }
  435.  
  436. /* ========================================================== */
  437. /* print matrix */
  438.  
  439. #define MAT_PRINT_3X3(mmm) {                \
  440.    int i,j;                        \
  441.    printf ("matrix mmm is \n");                \
  442.    if (mmm == NULL) {                    \
  443.       printf (" Null \n");                \
  444.    } else {                        \
  445.       for (i=0; i<3; i++) {                \
  446.          for (j=0; j<3; j++) {                \
  447.             printf ("%f ", mmm[i][j]);            \
  448.          }                        \
  449.          printf (" \n");                \
  450.       }                            \
  451.    }                            \
  452. }
  453.  
  454. /* ========================================================== */
  455. /* print matrix */
  456.  
  457. #define MAT_PRINT_2X3(mmm) {                \
  458.    int i,j;                        \
  459.    printf ("matrix mmm is \n");                \
  460.    if (mmm == NULL) {                    \
  461.       printf (" Null \n");                \
  462.    } else {                        \
  463.       for (i=0; i<2; i++) {                \
  464.          for (j=0; j<3; j++) {                \
  465.             printf ("%f ", mmm[i][j]);            \
  466.          }                        \
  467.          printf (" \n");                \
  468.       }                            \
  469.    }                            \
  470. }
  471.  
  472. /* ========================================================== */
  473. /* initialize matrix */
  474.  
  475. #define IDENTIFY_MATRIX_3X3(m)            \
  476. {                        \
  477.    m[0][0] = 1.0;                \
  478.    m[0][1] = 0.0;                \
  479.    m[0][2] = 0.0;                \
  480.                         \
  481.    m[1][0] = 0.0;                \
  482.    m[1][1] = 1.0;                \
  483.    m[1][2] = 0.0;                \
  484.                         \
  485.    m[2][0] = 0.0;                \
  486.    m[2][1] = 0.0;                \
  487.    m[2][2] = 1.0;                \
  488. }
  489.  
  490. /* ========================================================== */
  491. /* initialize matrix */
  492.  
  493. #define IDENTIFY_MATRIX_4X4(m)            \
  494. {                        \
  495.    m[0][0] = 1.0;                \
  496.    m[0][1] = 0.0;                \
  497.    m[0][2] = 0.0;                \
  498.    m[0][3] = 0.0;                \
  499.                         \
  500.    m[1][0] = 0.0;                \
  501.    m[1][1] = 1.0;                \
  502.    m[1][2] = 0.0;                \
  503.    m[1][3] = 0.0;                \
  504.                         \
  505.    m[2][0] = 0.0;                \
  506.    m[2][1] = 0.0;                \
  507.    m[2][2] = 1.0;                \
  508.    m[2][3] = 0.0;                \
  509.                         \
  510.    m[3][0] = 0.0;                \
  511.    m[3][1] = 0.0;                \
  512.    m[3][2] = 0.0;                \
  513.    m[3][3] = 1.0;                \
  514. }
  515.  
  516. /* ========================================================== */
  517. /* matrix copy */
  518.  
  519. #define COPY_MATRIX_2X2(b,a)    \
  520. {                \
  521.    b[0][0] = a[0][0];        \
  522.    b[0][1] = a[0][1];        \
  523.                 \
  524.    b[1][0] = a[1][0];        \
  525.    b[1][1] = a[1][1];        \
  526.                 \
  527. }
  528.  
  529. /* ========================================================== */
  530. /* matrix copy */
  531.  
  532. #define COPY_MATRIX_2X3(b,a)    \
  533. {                \
  534.    b[0][0] = a[0][0];        \
  535.    b[0][1] = a[0][1];        \
  536.    b[0][2] = a[0][2];        \
  537.                 \
  538.    b[1][0] = a[1][0];        \
  539.    b[1][1] = a[1][1];        \
  540.    b[1][2] = a[1][2];        \
  541. }
  542.  
  543. /* ========================================================== */
  544. /* matrix copy */
  545.  
  546. #define COPY_MATRIX_3X3(b,a)    \
  547. {                \
  548.    b[0][0] = a[0][0];        \
  549.    b[0][1] = a[0][1];        \
  550.    b[0][2] = a[0][2];        \
  551.                 \
  552.    b[1][0] = a[1][0];        \
  553.    b[1][1] = a[1][1];        \
  554.    b[1][2] = a[1][2];        \
  555.                 \
  556.    b[2][0] = a[2][0];        \
  557.    b[2][1] = a[2][1];        \
  558.    b[2][2] = a[2][2];        \
  559. }
  560.  
  561. /* ========================================================== */
  562. /* matrix copy */
  563.  
  564. #define COPY_MATRIX_4X4(b,a)    \
  565. {                \
  566.    b[0][0] = a[0][0];        \
  567.    b[0][1] = a[0][1];        \
  568.    b[0][2] = a[0][2];        \
  569.    b[0][3] = a[0][3];        \
  570.                 \
  571.    b[1][0] = a[1][0];        \
  572.    b[1][1] = a[1][1];        \
  573.    b[1][2] = a[1][2];        \
  574.    b[1][3] = a[1][3];        \
  575.                 \
  576.    b[2][0] = a[2][0];        \
  577.    b[2][1] = a[2][1];        \
  578.    b[2][2] = a[2][2];        \
  579.    b[2][3] = a[2][3];        \
  580.                 \
  581.    b[3][0] = a[3][0];        \
  582.    b[3][1] = a[3][1];        \
  583.    b[3][2] = a[3][2];        \
  584.    b[3][3] = a[3][3];        \
  585. }
  586.  
  587. /* ========================================================== */
  588. /* matrix transpose */
  589.  
  590. #define TRANSPOSE_MATRIX_2X2(b,a)    \
  591. {                \
  592.    b[0][0] = a[0][0];        \
  593.    b[0][1] = a[1][0];        \
  594.                 \
  595.    b[1][0] = a[0][1];        \
  596.    b[1][1] = a[1][1];        \
  597. }
  598.  
  599. /* ========================================================== */
  600. /* matrix transpose */
  601.  
  602. #define TRANSPOSE_MATRIX_3X3(b,a)    \
  603. {                \
  604.    b[0][0] = a[0][0];        \
  605.    b[0][1] = a[1][0];        \
  606.    b[0][2] = a[2][0];        \
  607.                 \
  608.    b[1][0] = a[0][1];        \
  609.    b[1][1] = a[1][1];        \
  610.    b[1][2] = a[2][1];        \
  611.                 \
  612.    b[2][0] = a[0][2];        \
  613.    b[2][1] = a[1][2];        \
  614.    b[2][2] = a[2][2];        \
  615. }
  616.  
  617. /* ========================================================== */
  618. /* matrix transpose */
  619.  
  620. #define TRANSPOSE_MATRIX_4X4(b,a)    \
  621. {                \
  622.    b[0][0] = a[0][0];        \
  623.    b[0][1] = a[1][0];        \
  624.    b[0][2] = a[2][0];        \
  625.    b[0][3] = a[3][0];        \
  626.                 \
  627.    b[1][0] = a[0][1];        \
  628.    b[1][1] = a[1][1];        \
  629.    b[1][2] = a[2][1];        \
  630.    b[1][3] = a[3][1];        \
  631.                 \
  632.    b[2][0] = a[0][2];        \
  633.    b[2][1] = a[1][2];        \
  634.    b[2][2] = a[2][2];        \
  635.    b[2][3] = a[3][2];        \
  636.                 \
  637.    b[3][0] = a[0][3];        \
  638.    b[3][1] = a[1][3];        \
  639.    b[3][2] = a[2][3];        \
  640.    b[3][3] = a[3][3];        \
  641. }
  642.  
  643. /* ========================================================== */
  644. /* multiply matrix by scalar */
  645.  
  646. #define SCALE_MATRIX_2X2(b,s,a)        \
  647. {                    \
  648.    b[0][0] = (s) * a[0][0];        \
  649.    b[0][1] = (s) * a[0][1];        \
  650.                     \
  651.    b[1][0] = (s) * a[1][0];        \
  652.    b[1][1] = (s) * a[1][1];        \
  653. }
  654.  
  655. /* ========================================================== */
  656. /* multiply matrix by scalar */
  657.  
  658. #define SCALE_MATRIX_3X3(b,s,a)        \
  659. {                    \
  660.    b[0][0] = (s) * a[0][0];        \
  661.    b[0][1] = (s) * a[0][1];        \
  662.    b[0][2] = (s) * a[0][2];        \
  663.                     \
  664.    b[1][0] = (s) * a[1][0];        \
  665.    b[1][1] = (s) * a[1][1];        \
  666.    b[1][2] = (s) * a[1][2];        \
  667.                     \
  668.    b[2][0] = (s) * a[2][0];        \
  669.    b[2][1] = (s) * a[2][1];        \
  670.    b[2][2] = (s) * a[2][2];        \
  671. }
  672.  
  673. /* ========================================================== */
  674. /* multiply matrix by scalar */
  675.  
  676. #define SCALE_MATRIX_4X4(b,s,a)        \
  677. {                    \
  678.    b[0][0] = (s) * a[0][0];        \
  679.    b[0][1] = (s) * a[0][1];        \
  680.    b[0][2] = (s) * a[0][2];        \
  681.    b[0][3] = (s) * a[0][3];        \
  682.                     \
  683.    b[1][0] = (s) * a[1][0];        \
  684.    b[1][1] = (s) * a[1][1];        \
  685.    b[1][2] = (s) * a[1][2];        \
  686.    b[1][3] = (s) * a[1][3];        \
  687.                     \
  688.    b[2][0] = (s) * a[2][0];        \
  689.    b[2][1] = (s) * a[2][1];        \
  690.    b[2][2] = (s) * a[2][2];        \
  691.    b[2][3] = (s) * a[2][3];        \
  692.                     \
  693.    b[3][0] = s * a[3][0];        \
  694.    b[3][1] = s * a[3][1];        \
  695.    b[3][2] = s * a[3][2];        \
  696.    b[3][3] = s * a[3][3];        \
  697. }
  698.  
  699. /* ========================================================== */
  700. /* multiply matrix by scalar */
  701.  
  702. #define ACCUM_SCALE_MATRIX_2X2(b,s,a)        \
  703. {                    \
  704.    b[0][0] += (s) * a[0][0];        \
  705.    b[0][1] += (s) * a[0][1];        \
  706.                     \
  707.    b[1][0] += (s) * a[1][0];        \
  708.    b[1][1] += (s) * a[1][1];        \
  709. }
  710.  
  711. /* +========================================================== */
  712. /* multiply matrix by scalar */
  713.  
  714. #define ACCUM_SCALE_MATRIX_3X3(b,s,a)        \
  715. {                    \
  716.    b[0][0] += (s) * a[0][0];        \
  717.    b[0][1] += (s) * a[0][1];        \
  718.    b[0][2] += (s) * a[0][2];        \
  719.                     \
  720.    b[1][0] += (s) * a[1][0];        \
  721.    b[1][1] += (s) * a[1][1];        \
  722.    b[1][2] += (s) * a[1][2];        \
  723.                     \
  724.    b[2][0] += (s) * a[2][0];        \
  725.    b[2][1] += (s) * a[2][1];        \
  726.    b[2][2] += (s) * a[2][2];        \
  727. }
  728.  
  729. /* +========================================================== */
  730. /* multiply matrix by scalar */
  731.  
  732. #define ACCUM_SCALE_MATRIX_4X4(b,s,a)        \
  733. {                    \
  734.    b[0][0] += (s) * a[0][0];        \
  735.    b[0][1] += (s) * a[0][1];        \
  736.    b[0][2] += (s) * a[0][2];        \
  737.    b[0][3] += (s) * a[0][3];        \
  738.                     \
  739.    b[1][0] += (s) * a[1][0];        \
  740.    b[1][1] += (s) * a[1][1];        \
  741.    b[1][2] += (s) * a[1][2];        \
  742.    b[1][3] += (s) * a[1][3];        \
  743.                     \
  744.    b[2][0] += (s) * a[2][0];        \
  745.    b[2][1] += (s) * a[2][1];        \
  746.    b[2][2] += (s) * a[2][2];        \
  747.    b[2][3] += (s) * a[2][3];        \
  748.                     \
  749.    b[3][0] += (s) * a[3][0];        \
  750.    b[3][1] += (s) * a[3][1];        \
  751.    b[3][2] += (s) * a[3][2];        \
  752.    b[3][3] += (s) * a[3][3];        \
  753. }
  754.  
  755. /* +========================================================== */
  756. /* matrix product */
  757. /* c[x][y] = a[x][0]*b[0][y]+a[x][1]*b[1][y]+a[x][2]*b[2][y]+a[x][3]*b[3][y];*/
  758.  
  759. #define MATRIX_PRODUCT_2X2(c,a,b)        \
  760. {                        \
  761.    c[0][0] = a[0][0]*b[0][0]+a[0][1]*b[1][0];    \
  762.    c[0][1] = a[0][0]*b[0][1]+a[0][1]*b[1][1];    \
  763.                         \
  764.    c[1][0] = a[1][0]*b[0][0]+a[1][1]*b[1][0];    \
  765.    c[1][1] = a[1][0]*b[0][1]+a[1][1]*b[1][1];    \
  766.                         \
  767. }
  768.  
  769. /* ========================================================== */
  770. /* matrix product */
  771. /* c[x][y] = a[x][0]*b[0][y]+a[x][1]*b[1][y]+a[x][2]*b[2][y]+a[x][3]*b[3][y];*/
  772.  
  773. #define MATRIX_PRODUCT_3X3(c,a,b)                \
  774. {                                \
  775.    c[0][0] = a[0][0]*b[0][0]+a[0][1]*b[1][0]+a[0][2]*b[2][0];    \
  776.    c[0][1] = a[0][0]*b[0][1]+a[0][1]*b[1][1]+a[0][2]*b[2][1];    \
  777.    c[0][2] = a[0][0]*b[0][2]+a[0][1]*b[1][2]+a[0][2]*b[2][2];    \
  778.                                 \
  779.    c[1][0] = a[1][0]*b[0][0]+a[1][1]*b[1][0]+a[1][2]*b[2][0];    \
  780.    c[1][1] = a[1][0]*b[0][1]+a[1][1]*b[1][1]+a[1][2]*b[2][1];    \
  781.    c[1][2] = a[1][0]*b[0][2]+a[1][1]*b[1][2]+a[1][2]*b[2][2];    \
  782.                                 \
  783.    c[2][0] = a[2][0]*b[0][0]+a[2][1]*b[1][0]+a[2][2]*b[2][0];    \
  784.    c[2][1] = a[2][0]*b[0][1]+a[2][1]*b[1][1]+a[2][2]*b[2][1];    \
  785.    c[2][2] = a[2][0]*b[0][2]+a[2][1]*b[1][2]+a[2][2]*b[2][2];    \
  786. }
  787.  
  788. /* ========================================================== */
  789. /* matrix product */
  790. /* c[x][y] = a[x][0]*b[0][y]+a[x][1]*b[1][y]+a[x][2]*b[2][y]+a[x][3]*b[3][y];*/
  791.  
  792. #define MATRIX_PRODUCT_4X4(c,a,b)        \
  793. {                        \
  794.    c[0][0] = a[0][0]*b[0][0]+a[0][1]*b[1][0]+a[0][2]*b[2][0]+a[0][3]*b[3][0];\
  795.    c[0][1] = a[0][0]*b[0][1]+a[0][1]*b[1][1]+a[0][2]*b[2][1]+a[0][3]*b[3][1];\
  796.    c[0][2] = a[0][0]*b[0][2]+a[0][1]*b[1][2]+a[0][2]*b[2][2]+a[0][3]*b[3][2];\
  797.    c[0][3] = a[0][0]*b[0][3]+a[0][1]*b[1][3]+a[0][2]*b[2][3]+a[0][3]*b[3][3];\
  798.                         \
  799.    c[1][0] = a[1][0]*b[0][0]+a[1][1]*b[1][0]+a[1][2]*b[2][0]+a[1][3]*b[3][0];\
  800.    c[1][1] = a[1][0]*b[0][1]+a[1][1]*b[1][1]+a[1][2]*b[2][1]+a[1][3]*b[3][1];\
  801.    c[1][2] = a[1][0]*b[0][2]+a[1][1]*b[1][2]+a[1][2]*b[2][2]+a[1][3]*b[3][2];\
  802.    c[1][3] = a[1][0]*b[0][3]+a[1][1]*b[1][3]+a[1][2]*b[2][3]+a[1][3]*b[3][3];\
  803.                         \
  804.    c[2][0] = a[2][0]*b[0][0]+a[2][1]*b[1][0]+a[2][2]*b[2][0]+a[2][3]*b[3][0];\
  805.    c[2][1] = a[2][0]*b[0][1]+a[2][1]*b[1][1]+a[2][2]*b[2][1]+a[2][3]*b[3][1];\
  806.    c[2][2] = a[2][0]*b[0][2]+a[2][1]*b[1][2]+a[2][2]*b[2][2]+a[2][3]*b[3][2];\
  807.    c[2][3] = a[2][0]*b[0][3]+a[2][1]*b[1][3]+a[2][2]*b[2][3]+a[2][3]*b[3][3];\
  808.                         \
  809.    c[3][0] = a[3][0]*b[0][0]+a[3][1]*b[1][0]+a[3][2]*b[2][0]+a[3][3]*b[3][0];\
  810.    c[3][1] = a[3][0]*b[0][1]+a[3][1]*b[1][1]+a[3][2]*b[2][1]+a[3][3]*b[3][1];\
  811.    c[3][2] = a[3][0]*b[0][2]+a[3][1]*b[1][2]+a[3][2]*b[2][2]+a[3][3]*b[3][2];\
  812.    c[3][3] = a[3][0]*b[0][3]+a[3][1]*b[1][3]+a[3][2]*b[2][3]+a[3][3]*b[3][3];\
  813. }
  814.  
  815. /* ========================================================== */
  816. /* matrix times vector */
  817.  
  818. #define MAT_DOT_VEC_2X2(p,m,v)                    \
  819. {                                \
  820.    p[0] = m[0][0]*v[0] + m[0][1]*v[1];                \
  821.    p[1] = m[1][0]*v[0] + m[1][1]*v[1];                \
  822. }
  823.  
  824. /* ========================================================== */
  825. /* matrix times vector */
  826.  
  827. #define MAT_DOT_VEC_3X3(p,m,v)                    \
  828. {                                \
  829.    p[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];        \
  830.    p[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];        \
  831.    p[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];        \
  832. }
  833.  
  834. /* ========================================================== */
  835. /* matrix times vector */
  836.  
  837. #define MAT_DOT_VEC_4X4(p,m,v)                    \
  838. {                                \
  839.    p[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2] + m[0][3]*v[3];    \
  840.    p[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2] + m[1][3]*v[3];    \
  841.    p[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2] + m[2][3]*v[3];    \
  842.    p[3] = m[3][0]*v[0] + m[3][1]*v[1] + m[3][2]*v[2] + m[3][3]*v[3];    \
  843. }
  844.  
  845. /* ========================================================== */
  846. /* vector transpose times matrix */
  847. /* p[j] = v[0]*m[0][j] + v[1]*m[1][j] + v[2]*m[2][j]; */
  848.  
  849. #define VEC_DOT_MAT_3X3(p,v,m)                    \
  850. {                                \
  851.    p[0] = v[0]*m[0][0] + v[1]*m[1][0] + v[2]*m[2][0];        \
  852.    p[1] = v[0]*m[0][1] + v[1]*m[1][1] + v[2]*m[2][1];        \
  853.    p[2] = v[0]*m[0][2] + v[1]*m[1][2] + v[2]*m[2][2];        \
  854. }
  855.  
  856. /* ========================================================== */
  857. /* affine matrix times vector */
  858. /* The matrix is assumed to be an affine matrix, with last two 
  859.  * entries representing a translation */
  860.  
  861. #define MAT_DOT_VEC_2X3(p,m,v)                    \
  862. {                                \
  863.    p[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2];        \
  864.    p[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2];        \
  865. }
  866.  
  867. /* ========================================================== */
  868. /* inverse transpose of matrix times vector
  869.  *
  870.  * This macro computes inverse transpose of matrix m, 
  871.  * and multiplies vector v into it, to yeild vector p
  872.  *
  873.  * DANGER !!! Do Not use this on normal vectors!!!
  874.  * It will leave normals the wrong length !!!
  875.  * See macro below for use on normals.
  876.  */
  877.  
  878. #define INV_TRANSP_MAT_DOT_VEC_2X2(p,m,v)            \
  879. {                                \
  880.    gleDouble det;                        \
  881.                                 \
  882.    det = m[0][0]*m[1][1] - m[0][1]*m[1][0];            \
  883.    p[0] = m[1][1]*v[0] - m[1][0]*v[1];                \
  884.    p[1] = - m[0][1]*v[0] + m[0][0]*v[1];            \
  885.                                 \
  886.    /* if matrix not singular, and not orthonormal, then renormalize */ \
  887.    if ((det!=1.0) && (det != 0.0)) {                \
  888.       det = 1.0 / det;                        \
  889.       p[0] *= det;                        \
  890.       p[1] *= det;                        \
  891.    }                                \
  892. }
  893.  
  894. /* ========================================================== */
  895. /* transform normal vector by inverse transpose of matrix 
  896.  * and then renormalize the vector 
  897.  *
  898.  * This macro computes inverse transpose of matrix m, 
  899.  * and multiplies vector v into it, to yeild vector p
  900.  * Vector p is then normalized.
  901.  */
  902.  
  903.  
  904. #define NORM_XFORM_2X2(p,m,v)                    \
  905. {                                \
  906.    double len;                            \
  907.                                 \
  908.    /* do nothing if off-diagonals are zero and diagonals are     \
  909.     * equal */                            \
  910.    if ((m[0][1] != 0.0) || (m[1][0] != 0.0) || (m[0][0] != m[1][1])) { \
  911.       p[0] = m[1][1]*v[0] - m[1][0]*v[1];            \
  912.       p[1] = - m[0][1]*v[0] + m[0][0]*v[1];            \
  913.                                 \
  914.       len = p[0]*p[0] + p[1]*p[1];                \
  915.       len = 1.0 / sqrt (len);                    \
  916.       p[0] *= len;                        \
  917.       p[1] *= len;                        \
  918.    } else {                            \
  919.       VEC_COPY_2 (p, v);                    \
  920.    }                                \
  921. }
  922.  
  923. /* ========================================================== */
  924. /* outer product of vector times vector transpose 
  925.  *
  926.  * The outer product of vector v and vector transpose t yeilds 
  927.  * dyadic matrix m.
  928.  */
  929.  
  930. #define OUTER_PRODUCT_2X2(m,v,t)                \
  931. {                                \
  932.    m[0][0] = v[0] * t[0];                    \
  933.    m[0][1] = v[0] * t[1];                    \
  934.                                 \
  935.    m[1][0] = v[1] * t[0];                    \
  936.    m[1][1] = v[1] * t[1];                    \
  937. }
  938.  
  939. /* ========================================================== */
  940. /* outer product of vector times vector transpose 
  941.  *
  942.  * The outer product of vector v and vector transpose t yeilds 
  943.  * dyadic matrix m.
  944.  */
  945.  
  946. #define OUTER_PRODUCT_3X3(m,v,t)                \
  947. {                                \
  948.    m[0][0] = v[0] * t[0];                    \
  949.    m[0][1] = v[0] * t[1];                    \
  950.    m[0][2] = v[0] * t[2];                    \
  951.                                 \
  952.    m[1][0] = v[1] * t[0];                    \
  953.    m[1][1] = v[1] * t[1];                    \
  954.    m[1][2] = v[1] * t[2];                    \
  955.                                 \
  956.    m[2][0] = v[2] * t[0];                    \
  957.    m[2][1] = v[2] * t[1];                    \
  958.    m[2][2] = v[2] * t[2];                    \
  959. }
  960.  
  961. /* ========================================================== */
  962. /* outer product of vector times vector transpose 
  963.  *
  964.  * The outer product of vector v and vector transpose t yeilds 
  965.  * dyadic matrix m.
  966.  */
  967.  
  968. #define OUTER_PRODUCT_4X4(m,v,t)                \
  969. {                                \
  970.    m[0][0] = v[0] * t[0];                    \
  971.    m[0][1] = v[0] * t[1];                    \
  972.    m[0][2] = v[0] * t[2];                    \
  973.    m[0][3] = v[0] * t[3];                    \
  974.                                 \
  975.    m[1][0] = v[1] * t[0];                    \
  976.    m[1][1] = v[1] * t[1];                    \
  977.    m[1][2] = v[1] * t[2];                    \
  978.    m[1][3] = v[1] * t[3];                    \
  979.                                 \
  980.    m[2][0] = v[2] * t[0];                    \
  981.    m[2][1] = v[2] * t[1];                    \
  982.    m[2][2] = v[2] * t[2];                    \
  983.    m[2][3] = v[2] * t[3];                    \
  984.                                 \
  985.    m[3][0] = v[3] * t[0];                    \
  986.    m[3][1] = v[3] * t[1];                    \
  987.    m[3][2] = v[3] * t[2];                    \
  988.    m[3][3] = v[3] * t[3];                    \
  989. }
  990.  
  991. /* +========================================================== */
  992. /* outer product of vector times vector transpose 
  993.  *
  994.  * The outer product of vector v and vector transpose t yeilds 
  995.  * dyadic matrix m.
  996.  */
  997.  
  998. #define ACCUM_OUTER_PRODUCT_2X2(m,v,t)                \
  999. {                                \
  1000.    m[0][0] += v[0] * t[0];                    \
  1001.    m[0][1] += v[0] * t[1];                    \
  1002.                                 \
  1003.    m[1][0] += v[1] * t[0];                    \
  1004.    m[1][1] += v[1] * t[1];                    \
  1005. }
  1006.  
  1007. /* +========================================================== */
  1008. /* outer product of vector times vector transpose 
  1009.  *
  1010.  * The outer product of vector v and vector transpose t yeilds 
  1011.  * dyadic matrix m.
  1012.  */
  1013.  
  1014. #define ACCUM_OUTER_PRODUCT_3X3(m,v,t)                \
  1015. {                                \
  1016.    m[0][0] += v[0] * t[0];                    \
  1017.    m[0][1] += v[0] * t[1];                    \
  1018.    m[0][2] += v[0] * t[2];                    \
  1019.                                 \
  1020.    m[1][0] += v[1] * t[0];                    \
  1021.    m[1][1] += v[1] * t[1];                    \
  1022.    m[1][2] += v[1] * t[2];                    \
  1023.                                 \
  1024.    m[2][0] += v[2] * t[0];                    \
  1025.    m[2][1] += v[2] * t[1];                    \
  1026.    m[2][2] += v[2] * t[2];                    \
  1027. }
  1028.  
  1029. /* +========================================================== */
  1030. /* outer product of vector times vector transpose 
  1031.  *
  1032.  * The outer product of vector v and vector transpose t yeilds 
  1033.  * dyadic matrix m.
  1034.  */
  1035.  
  1036. #define ACCUM_OUTER_PRODUCT_4X4(m,v,t)                \
  1037. {                                \
  1038.    m[0][0] += v[0] * t[0];                    \
  1039.    m[0][1] += v[0] * t[1];                    \
  1040.    m[0][2] += v[0] * t[2];                    \
  1041.    m[0][3] += v[0] * t[3];                    \
  1042.                                 \
  1043.    m[1][0] += v[1] * t[0];                    \
  1044.    m[1][1] += v[1] * t[1];                    \
  1045.    m[1][2] += v[1] * t[2];                    \
  1046.    m[1][3] += v[1] * t[3];                    \
  1047.                                 \
  1048.    m[2][0] += v[2] * t[0];                    \
  1049.    m[2][1] += v[2] * t[1];                    \
  1050.    m[2][2] += v[2] * t[2];                    \
  1051.    m[2][3] += v[2] * t[3];                    \
  1052.                                 \
  1053.    m[3][0] += v[3] * t[0];                    \
  1054.    m[3][1] += v[3] * t[1];                    \
  1055.    m[3][2] += v[3] * t[2];                    \
  1056.    m[3][3] += v[3] * t[3];                    \
  1057. }
  1058.  
  1059. /* +========================================================== */
  1060. /* determinant of matrix
  1061.  *
  1062.  * Computes determinant of matrix m, returning d
  1063.  */
  1064.  
  1065. #define DETERMINANT_2X2(d,m)                    \
  1066. {                                \
  1067.    d = m[0][0] * m[1][1] - m[0][1] * m[1][0];            \
  1068. }
  1069.  
  1070. /* ========================================================== */
  1071. /* determinant of matrix
  1072.  *
  1073.  * Computes determinant of matrix m, returning d
  1074.  */
  1075.  
  1076. #define DETERMINANT_3X3(d,m)                    \
  1077. {                                \
  1078.    d = m[0][0] * (m[1][1]*m[2][2] - m[1][2] * m[2][1]);        \
  1079.    d -= m[0][1] * (m[1][0]*m[2][2] - m[1][2] * m[2][0]);    \
  1080.    d += m[0][2] * (m[1][0]*m[2][1] - m[1][1] * m[2][0]);    \
  1081. }
  1082.  
  1083. /* ========================================================== */
  1084. /* i,j,th cofactor of a 4x4 matrix
  1085.  *
  1086.  */
  1087.  
  1088. #define COFACTOR_4X4_IJ(fac,m,i,j)                 \
  1089. {                                \
  1090.    int ii[4], jj[4], k;                        \
  1091.                                 \
  1092.    /* compute which row, columnt to skip */            \
  1093.    for (k=0; k<i; k++) ii[k] = k;                \
  1094.    for (k=i; k<3; k++) ii[k] = k+1;                \
  1095.    for (k=0; k<j; k++) jj[k] = k;                \
  1096.    for (k=j; k<3; k++) jj[k] = k+1;                \
  1097.                                 \
  1098.    (fac) = m[ii[0]][jj[0]] * (m[ii[1]][jj[1]]*m[ii[2]][jj[2]]     \
  1099.                             - m[ii[1]][jj[2]]*m[ii[2]][jj[1]]); \
  1100.    (fac) -= m[ii[0]][jj[1]] * (m[ii[1]][jj[0]]*m[ii[2]][jj[2]]    \
  1101.                              - m[ii[1]][jj[2]]*m[ii[2]][jj[0]]);\
  1102.    (fac) += m[ii[0]][jj[2]] * (m[ii[1]][jj[0]]*m[ii[2]][jj[1]]    \
  1103.                              - m[ii[1]][jj[1]]*m[ii[2]][jj[0]]);\
  1104.                                 \
  1105.    /* compute sign */                        \
  1106.    k = i+j;                            \
  1107.    if ( k != (k/2)*2) {                        \
  1108.       (fac) = -(fac);                        \
  1109.    }                                \
  1110. }
  1111.  
  1112. /* ========================================================== */
  1113. /* determinant of matrix
  1114.  *
  1115.  * Computes determinant of matrix m, returning d
  1116.  */
  1117.  
  1118. #define DETERMINANT_4X4(d,m)                    \
  1119. {                                \
  1120.    double cofac;                        \
  1121.    COFACTOR_4X4_IJ (cofac, m, 0, 0);                \
  1122.    d = m[0][0] * cofac;                        \
  1123.    COFACTOR_4X4_IJ (cofac, m, 0, 1);                \
  1124.    d += m[0][1] * cofac;                    \
  1125.    COFACTOR_4X4_IJ (cofac, m, 0, 2);                \
  1126.    d += m[0][2] * cofac;                    \
  1127.    COFACTOR_4X4_IJ (cofac, m, 0, 3);                \
  1128.    d += m[0][3] * cofac;                    \
  1129. }
  1130.  
  1131. /* ========================================================== */
  1132. /* cofactor of matrix
  1133.  *
  1134.  * Computes cofactor of matrix m, returning a
  1135.  */
  1136.  
  1137. #define COFACTOR_2X2(a,m)                    \
  1138. {                                \
  1139.    a[0][0] = (m)[1][1];                        \
  1140.    a[0][1] = - (m)[1][0];                        \
  1141.    a[1][0] = - (m)[0][1];                        \
  1142.    a[1][1] = (m)[0][0];                        \
  1143. }
  1144.  
  1145. /* ========================================================== */
  1146. /* cofactor of matrix
  1147.  *
  1148.  * Computes cofactor of matrix m, returning a
  1149.  */
  1150.  
  1151. #define COFACTOR_3X3(a,m)                    \
  1152. {                                \
  1153.    a[0][0] = m[1][1]*m[2][2] - m[1][2]*m[2][1];            \
  1154.    a[0][1] = - (m[1][0]*m[2][2] - m[2][0]*m[1][2]);        \
  1155.    a[0][2] = m[1][0]*m[2][1] - m[1][1]*m[2][0];            \
  1156.    a[1][0] = - (m[0][1]*m[2][2] - m[0][2]*m[2][1]);        \
  1157.    a[1][1] = m[0][0]*m[2][2] - m[0][2]*m[2][0];            \
  1158.    a[1][2] = - (m[0][0]*m[2][1] - m[0][1]*m[2][0]);        \
  1159.    a[2][0] = m[0][1]*m[1][2] - m[0][2]*m[1][1];            \
  1160.    a[2][1] = - (m[0][0]*m[1][2] - m[0][2]*m[1][0]);        \
  1161.    a[2][2] = m[0][0]*m[1][1] - m[0][1]*m[1][0]);        \
  1162. }
  1163.  
  1164. /* ========================================================== */
  1165. /* cofactor of matrix
  1166.  *
  1167.  * Computes cofactor of matrix m, returning a
  1168.  */
  1169.  
  1170. #define COFACTOR_4X4(a,m)                    \
  1171. {                                \
  1172.    int i,j;                            \
  1173.                                 \
  1174.    for (i=0; i<4; i++) {                    \
  1175.       for (j=0; j<4; j++) {                    \
  1176.          COFACTOR_4X4_IJ (a[i][j], m, i, j);            \
  1177.       }                                \
  1178.    }                                \
  1179. }
  1180.  
  1181. /* ========================================================== */
  1182. /* adjoint of matrix
  1183.  *
  1184.  * Computes adjoint of matrix m, returning a
  1185.  * (Note that adjoint is just the transpose of the cofactor matrix)
  1186.  */
  1187.  
  1188. #define ADJOINT_2X2(a,m)                    \
  1189. {                                \
  1190.    a[0][0] = (m)[1][1];                        \
  1191.    a[1][0] = - (m)[1][0];                        \
  1192.    a[0][1] = - (m)[0][1];                        \
  1193.    a[1][1] = (m)[0][0];                        \
  1194. }
  1195.  
  1196. /* ========================================================== */
  1197. /* adjoint of matrix
  1198.  *
  1199.  * Computes adjoint of matrix m, returning a
  1200.  * (Note that adjoint is just the transpose of the cofactor matrix)
  1201.  */
  1202.  
  1203.  
  1204. #define ADJOINT_3X3(a,m)                    \
  1205. {                                \
  1206.    a[0][0] = m[1][1]*m[2][2] - m[1][2]*m[2][1];            \
  1207.    a[1][0] = - (m[1][0]*m[2][2] - m[2][0]*m[1][2]);        \
  1208.    a[2][0] = m[1][0]*m[2][1] - m[1][1]*m[2][0];            \
  1209.    a[0][1] = - (m[0][1]*m[2][2] - m[0][2]*m[2][1]);        \
  1210.    a[1][1] = m[0][0]*m[2][2] - m[0][2]*m[2][0];            \
  1211.    a[2][1] = - (m[0][0]*m[2][1] - m[0][1]*m[2][0]);        \
  1212.    a[0][2] = m[0][1]*m[1][2] - m[0][2]*m[1][1];            \
  1213.    a[1][2] = - (m[0][0]*m[1][2] - m[0][2]*m[1][0]);        \
  1214.    a[2][2] = m[0][0]*m[1][1] - m[0][1]*m[1][0]);        \
  1215. }
  1216.  
  1217. /* ========================================================== */
  1218. /* adjoint of matrix
  1219.  *
  1220.  * Computes adjoint of matrix m, returning a
  1221.  * (Note that adjoint is just the transpose of the cofactor matrix)
  1222.  */
  1223.  
  1224. #define ADJOINT_4X4(a,m)                    \
  1225. {                                \
  1226.    int i,j;                            \
  1227.                                 \
  1228.    for (i=0; i<4; i++) {                    \
  1229.       for (j=0; j<4; j++) {                    \
  1230.          COFACTOR_4X4_IJ (a[j][i], m, i, j);            \
  1231.       }                                \
  1232.    }                                \
  1233. }
  1234.  
  1235. /* ========================================================== */
  1236. /* compute adjoint of matrix and scale
  1237.  *
  1238.  * Computes adjoint of matrix m, scales it by s, returning a
  1239.  */
  1240.  
  1241. #define SCALE_ADJOINT_2X2(a,s,m)                \
  1242. {                                \
  1243.    a[0][0] = (s) * m[1][1];                    \
  1244.    a[1][0] = - (s) * m[1][0];                    \
  1245.    a[0][1] = - (s) * m[0][1];                    \
  1246.    a[1][1] = (s) * m[0][0];                    \
  1247. }
  1248.  
  1249. /* ========================================================== */
  1250. /* compute adjoint of matrix and scale
  1251.  *
  1252.  * Computes adjoint of matrix m, scales it by s, returning a
  1253.  */
  1254.  
  1255. #define SCALE_ADJOINT_3X3(a,s,m)                \
  1256. {                                \
  1257.    a[0][0] = (s) * (m[1][1] * m[2][2] - m[1][2] * m[2][1]);    \
  1258.    a[1][0] = (s) * (m[1][2] * m[2][0] - m[1][0] * m[2][2]);    \
  1259.    a[2][0] = (s) * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);    \
  1260.                                 \
  1261.    a[0][1] = (s) * (m[0][2] * m[2][1] - m[0][1] * m[2][2]);    \
  1262.    a[1][1] = (s) * (m[0][0] * m[2][2] - m[0][2] * m[2][0]);    \
  1263.    a[2][1] = (s) * (m[0][1] * m[2][0] - m[0][0] * m[2][1]);    \
  1264.                                 \
  1265.    a[0][2] = (s) * (m[0][1] * m[1][2] - m[0][2] * m[1][1]);    \
  1266.    a[1][2] = (s) * (m[0][2] * m[1][0] - m[0][0] * m[1][2]);    \
  1267.    a[2][2] = (s) * (m[0][0] * m[1][1] - m[0][1] * m[1][0]);    \
  1268. }
  1269.  
  1270. /* ========================================================== */
  1271. /* compute adjoint of matrix and scale
  1272.  *
  1273.  * Computes adjoint of matrix m, scales it by s, returning a
  1274.  */
  1275.  
  1276. #define SCALE_ADJOINT_4X4(a,s,m)                \
  1277. {                                \
  1278.    int i,j;                            \
  1279.                                 \
  1280.    for (i=0; i<4; i++) {                    \
  1281.       for (j=0; j<4; j++) {                    \
  1282.          COFACTOR_4X4_IJ (a[j][i], m, i, j);            \
  1283.          a[j][i] *= s;                        \
  1284.       }                                \
  1285.    }                                \
  1286. }
  1287.  
  1288. /* ========================================================== */
  1289. /* inverse of matrix 
  1290.  *
  1291.  * Compute inverse of matrix a, returning determinant m and 
  1292.  * inverse b
  1293.  */
  1294.  
  1295. #define INVERT_2X2(b,det,a)            \
  1296. {                        \
  1297.    double tmp;                    \
  1298.    DETERMINANT_2X2 (det, a);            \
  1299.    tmp = 1.0 / (det);                \
  1300.    SCALE_ADJOINT_2X2 (b, tmp, a);        \
  1301. }
  1302.  
  1303. /* ========================================================== */
  1304. /* inverse of matrix 
  1305.  *
  1306.  * Compute inverse of matrix a, returning determinant m and 
  1307.  * inverse b
  1308.  */
  1309.  
  1310. #define INVERT_3X3(b,det,a)            \
  1311. {                        \
  1312.    double tmp;                    \
  1313.    DETERMINANT_3X3 (det, a);            \
  1314.    tmp = 1.0 / (det);                \
  1315.    SCALE_ADJOINT_3X3 (b, tmp, a);        \
  1316. }
  1317.  
  1318. /* ========================================================== */
  1319. /* inverse of matrix 
  1320.  *
  1321.  * Compute inverse of matrix a, returning determinant m and 
  1322.  * inverse b
  1323.  */
  1324.  
  1325. #define INVERT_4X4(b,det,a)            \
  1326. {                        \
  1327.    double tmp;                    \
  1328.    DETERMINANT_4X4 (det, a);            \
  1329.    tmp = 1.0 / (det);                \
  1330.    SCALE_ADJOINT_4X4 (b, tmp, a);        \
  1331. }
  1332.  
  1333. /* ========================================================== */
  1334. #if defined(__cplusplus) || defined(c_plusplus)
  1335. }
  1336. #endif
  1337.  
  1338. #endif /* __GUTIL_VECTOR_H__ */
  1339. /* ===================== END OF FILE ======================== */
  1340.